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Abstract 

The correction of the self-interaction error (SIE) that is inherent to all standard density functional 
theory (DFT) calculations is an object of increasing interest. In this article we apply the very 
recently developed Fermi-orbital based approach for the self-interaction correction (FOSIC) [U |2] 
to a set of different molecular systems. Our study covers systems ranging from simple diatomic to 
large organic molecules. We focus our analysis on the direct estimation of the ionization potential 
from orbital eigenvalues. 
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I. INTRODUCTION 


Electronic structure calculations within the DFT framework have become an invaluable 
tool for physicists, chemists and material scientists due to low computational costs and 
sufficient accuracy of this method j3]. Within DFT, the quafity of the resufts and therefore 
of subsequent derived properties is strictfy bound to the functionaf used to evaluate the 
exchange and correlation energy (Exc)- The estimation of Exc is the only approximation 
that is needed in the expression of the total energy within DFT [Ij. 

The search for improved descriptions of this contribution has lead from the simple LDA 
and GGA approaches to various levels of hybrid functionals for example B3LYP, HSE and 
meta hybrid GGA’s like M06 [SHTj. 

Some of the failures of the different functionals have been considered to be intrinsic to the 
DFT approach. Examples are the dissociation energies of two center three electron systems 
[H], overestimation of the magnetic coupling P-TiTj. description of charge-transfer systems 
[niiiH] and ionization potentials OHS]. In nearly all cases these faults can be related to 
the self-interaction error which is, of course to a different degree, inherent in all functionals 
used in todays DFT implementations. This effect arises from the spurious interaction of an 
electron with itself. In the Hartee-Fock framework this contribution is totally cancelled by 
the exchange contribution and this is also the main reason why the hybrid functionals often 
predict properties closer to experiment [Hj. The self-interaction error is partially cancelled 
by the fraction of HF density that is used in the construction of Exc- 

The self-interaction problem within the DFT approach has been identified very early, 
already in the process the method was developed. Perdew and Zunger (PZ) proposed a 
method for self-interaction correction (SIC) [I7| back in 1981. The major drawback of 
the so called PZ-SIC method is, that one needs to evaluate an orbital dependent exchange 
correlation potential. This results in excessive computational cost and therefore in practice 
PZ-SIC is only applied to small model systems or simplified Hamiltonians [THH2n|. Some 
other approaches were developed and the focus was always to somehow avoid the calculation 
of the orbital dependent Exc [HI 122] • 

The aim of this paper is to investigate the performance of the FOSIC approach in the field 
of molecular systems. We focus on the evaluation of the vertical ionization potentials (IPs) 
for test cases, ranging from very small to large conjugated molecular systems. Although the 
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first IPs could be directly calculated using the DFT extension of the Koopmans’ theorem 
PSl El] the negative of the HOMO energy is typically several electron volts too small with 
respect to the experimental values. This effect, due to orbital self-interaction error, varies 
with the physical extent of the occupied and unoccupied orbitals and can be responsible for 
non-systematic errors that cause unphysical charge-transfer between chemically disparate 
components in molecular assemblies and solids. With the inclusion of the self-interaction 
correction the HOMO energies become much closer to the experimental IPs. Due to the fact 
that this property is directly related to the self-interaction error this provides a benchmark 
which is essential to validate the FOSIC method. It should be noted here that the Koop¬ 
mans’ theorem in Hartree Fock (HF) states that the HF eigenvalue is equal to the difference 
between the fully relaxed ground-state of the N electron system and the ground-state of the 
N-1 electron system when the energy of the N-1 electron system is relaxed within the space 
spanned by the orbitals of the N-electron system. Hence in HF, where the Koopmans’ 
theorem is exact, one should expect that the HF eigenvalue overestimates the fully 
relaxed HF A-SCF ionization energy and, in the weak-correlation limit, the experimental 
ionization energy as well. A similar variational Koopmans’ theorem has been proven for 
the PZ-SIC theory m but it has been noted that the non-quadratic dependence of the 
exchange-correlation functional removes the algebraically exact differences between the total 
energy of the N and N-1 electron and requires small, but difficult to account for, deviations 
between the eigenvalue and the energy difference. These differences have been referred to 
in the past [2S| as ” non-Koopmans” corrections and again more recently by Dabo et al 
[26l [27] . In contrast to HF, even in the weak-correlation limit, the analytically dependence 
of the energy functional on density prevents one from simply stating that the SIC-DFT 
eigenvalue should always overestimate the experimental ionization potential. However since 
to lowest order the SIC-DFT eigenvalue does not include relaxation of the system out of the 
space of orbitals spanned by the N-electron ground state, it may be reasonable to antici¬ 
pate that SIC-DFT eigenvalues which overestimate experiment are more likely in most cases. 

The paper is organized as follows: We will first give a short review on the essential ideas 
of SIC and the FOSIC approach to DFT. We will only focus on the general concepts that 
are necessary for this paper. For a detailed review we point the reader to the original papers 
published by Pederson and coworkers and the references given therein Second, we will 
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describe the details of the used calculation procedure. Finally we will discuss the obtained 
results and identify crucial points for the further development of the method. 


II. THEORETICAL BACKGROUND 


A. Fermi Orbital SIC 


The FOSIC method is based on the original Perdew and Zunger approach [28] . The DFT 
functional is expressed in the following way: 

nt] - (1) 


with 


J^PZ-SIC = _ {[/ J „ 0|} , (2) 

In equationthe orbitals {(paa} define orbital densities according to na<T(r) ~ 

The term U [pa,a] is the exact self-Coulomb energy and [Pa,cr,0] is the approximated 

exchange-correlation energy. Essentially the PZ-SIC method is a correction to the standard 
Exc given by the local spin density approximation. As discussed in [SHHEI] the best option 
to solve equation is the use of localized, canonical orbitals. Thereby one has to make 
sure that the Hermitian Lagrange multiplier matrix of the orbitals has to fullfill a O(A^) 
localization equation 






SIC 


(3) 


with respect to the orbital density Uaa- 

In principle any means for parameterizing a unitary transformation may be used to solve 
the localization equations. However, it has been recognized that the full optimization of the 
unitary transformation, at least when used within existing functionals, leads to an expression 
for the total energy that is not necessarily size extensive . Ensuring size extensivity within the 
PZ-SIC method is difficult since it is hard to systematically derive approximate functionals 
that will deliver a negative SIC energy for the highest occupied orbitals of a heavy atom. 
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A new approach, which constrains the functional to be unitarily invariant and fixes the 
problem associated with size extensivity, was to use localized orbitals on the basis of the 
Fermi hole the so called Fermi orbitals [S2l ES] • 

A set of KS orbitals can be transformed to Fermi orbitals in any point of space according 
to 


F^a{v) = 
FUv) = 


\/Per (^icr) 

'0acr(^icr)'0acr(r) 




(4) 

(5) 


Due to the procedure proposed in [T] it is possible to bypass the direct solution of the 
localization equations and therefore the computational effort needed has the advantage 
to have the scaling of standard DFT calculations. The expression for the exact exchange 
energy using a single Slater determinant is given by 

Ea 


= J dh j d^J 


( 6 ) 


According to P this can be rewritten to 




1 
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p^(r,a)p^(a, r) 
\/P(T(r)y/Pa(r) 



(7) 


where the term 

p^(r,a)p^(a,r) ^ .p. xp 


( 8 ) 


is the exchange hole and therefore equivalent to the square of the Fermi orbital. The con¬ 
structed Fermi orbitals depend parametrically on a set of N^- positions which correspond 
to qnasi-classical electronic positions. By the availability of gradients of the Fermi orbitals 
[2] one has a way to minimize the self interaction corrected total energy as function of the 
positions of the Fermi orbitals. 


B. Obtaining the Fermi orbital positions 

An important challenge of the FOSIC method is that one needs a set of initial Fermi 
orbitals in the beginning of the calculation. For atoms and very small molecules this could 
be obtained by a brute-force method (see [1]) or by using Themical intuition’. Clearly this 
may become not feasible for large molecules or automated calculations. 
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As a straightforward way to obtain approximate initial positions for the Fermi orbitals 
we propose using maximally localized Wannier orbitals (MLWF) |Hlj. It has been claimed 
that Wannier orbitals and Fermi orbitals are identical under certain conditions |SS] and they 
seem therefore as a natural choice to develop an initial-guess. Here we followed the approach 
of [36] and considered a set of Nocc eigenstates of The total energy is invariant with 

respect to unitary transformations Umn among the eigenstates 

NOCC 

|n^n) ^ ^ Umn l^m) • (9) 

m=l 

The unitary matrix U is chosen such that the resulting Ngcc orbitals {wn{r)} minimize their 
total quadratic spread, given by 

^n = Yl {{'^n\r‘^\Wn) “ (Wn|r|Wn)^) = ’ (^O) 

n n 

Each MLWF is characterized by a value of its quadratic spread, and its centre r„. 

We carried out the following procedure to obtain the initial guess for the Fermi orbital 
positions: 

1. perform a standard DFT calculation of the molecule of interest 

2. construct a set of MLWF’s according to equations and 

3. take the centers as initial positions for the subsequent FOSIC calculation 

These initial positions then need to be optimized in order to minimize Etot of the system 
of interest. In our investigation it turned out, that the is a quite robust initial guess. 
In Fig. [^) and b) we show for one example the initial guess and the final Fermi orbital 
positions after the optimization procedure. The quality of the guess is in general reasonable 
and reflects the qualitative details of the final Fermi orbital positions. However, we observed 
that with the current localization scheme the resulting Wannier centers appear to be too 
close to the nearest neighbour atoms. This is not a principle problem since the subsequent 
optimization of these initial positions always leads to correct placement of the final Fermi 
orbitals (see Fig. [^. A drawback is that in the current implementation the optimization 
algorithm needs a large number of steps to find the final positions of the Fermi orbitals. 
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FIG. 1: Example of the optimization process of the FO positions: a) The red markers show the 
initial positions of the FO centroids as obtained from the Wannier analysis, b) The green markers 
show the final converged positions of the FO centroids, c) Development of the total energy and 
the forces on the FO’s as a function of optimization step. A version of the FIRE algorithm m 
was used. 

C. Computational details 

We have used a modified version of the nrlmol all-electron DFT code [3HI EH] that uses 
large Gaussian-orbital basis sets |10] for the representation of the electron wavefunctions. 
The standard DFT potentials as well as the FOSIC potential is calculated on a mesh using 
the variational mesh technique [H]. The PW92 functional [l2] within the local density 
approximation (LDA) was used for the standard LDA as well as for the FOSIC calculations. 

The GGA calculations utilizing the B3LYP functional for exchange and correlation m 
sg were carried out using version 2.9 of the ORCA DFT package [IS]. We utilized the 
def2-TZVP basis set [IHI in order to avoid artifacts related to small basis sets. 

The geometry of all molecules have been optimized either with LDA/PW92 or the 
GGA/B3LYP functional to obtain relaxed structures with atomic forces below 0.01 eV/A. 

D. Ionization potential 

The primary ionization potentials for LDA and B3LYP were taken as the negative of 
the eigenvalue of the highest occupied molecular orbital IP = —chomo (see Table |T] and B- 
We further used the LDA/PW92 result as our starting point for the FOSIC calculations. 
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Following the description in sec. |IIB| we constructed a set of maximally localized Wannier 
orbitals from the Kohn-Sham orbitals of the LDA/PW92 calculation. The centers of the 
Wannier orbitals were used as initial positions for Fermi orbitals. A subsequent optimization 
of the Fermi orbital positions in order to minimize the total energy was carried out by using a 
conjugated gradient method. For this task the forces on the Fermi orbitals were calculated. 
The optimization algorithm was stopped once these forces on electrons were below 10“^ 
Hartree/bohr and —chomo was taken from the resulting electronic configuration. 

The ASCF values are obtained by taking the difference of the total energies between the 
neutral and the positive charged molecules without relaxing the geometry of the charged 
state. 


III. RESULTS AND DISCUSSION 

We have chosen two different groups of molecules in order to test how the FOSIC 
method performs within structures of different complexity. Test set 1 (see table |T| includes 
very small di- and triatomic molecules and small to medium sized molecules with double 
and triple bonds covering most of the main group elements. The 2nd test set consists of 
different organic aromatic structures ranging from simple 5 or 6 member rings with different 
substitutes to conjugated structures with increasing size (from two up to five condensed 
rings). 

As expected the uncorrected —ch values (no FOSIC) for the test set of very small 
molecules have a large mean average error (MAE) of 4.6 eV. The experimental ioniza¬ 
tion potential is generally underestimated my several electron volts. A much smaller error 
of 0.4 eV is found for the ASCF approach. The values do agree very well with the ex¬ 
perimental data due to the fact that for ASCF the self-interaction error cancels at least 
partially. In some examples like CO the ASCF-IP matches the experimental value almost 
perfectly. The —ch values obtained by the FOSIC approach on the other hand offer IP 
values with a mean average error of 2.2 eV. That is a clear improvement compared to the 
uncorrected LDA values. The — ch values obtained for the popular B3LYP functional have a 
mean average error of 3.6 eV which is significant larger than the FOSIC error. The B3LYP 
IP values are in between the LDA and FOSIC values. Both LDA and B3LYP underestimate 



the experimental ionization potential values while FOSIC typically overestimates. 

For the medium to large molecules presented in table [TT] no clear trend is visible. Overall 
the obtained —ch values for the FOSIC calculations have a lower mean average error of 1.4 eV 
compared to 1.9 eV the standard DFT values. For B3LYP (MAE = 2.3 eV) one notes again 
that FOSIC does perform better in general. Again we see that FOSIC overestimates the 
ionization potential. However the magnitude of the overestimation is different. The deviation 
of the FOSIC IP’s from the experimental values does not show a systematic behavior. All the 
differently calculated IP values, in comparison to there experimental counter part are shown 
in figure However, it is worth noting that especially for the larger molecules the deviation 
from experiment is rather small compared to LDA or B3LYP. in the case of chrysene and 
picene FOSIC provides even a qualitative better result. While LDA gives nearly identical 
IP’s for the two molecules FOSIC does reproduce the experimental finding, that picene has 

TABLE I: Test set 1: Experimental values for the experimental ionization potential (IP) and the 
respective calculated negative of the HOMO energy for different methods (all values in [eV]). MAE 
is the mean average error. As discussed in the text, an exact Koopmans’ theorem should overesti¬ 
mate the ionization energy. It is only exact in the limit of weak correlation, a quadratic expression 
for the exchange-correlation energy, and if electronic relaxation, outside the space spanned by the 
N-electron space of orbitals, is unimportant. 


molecule 

exp. IP 

^LDA 

^FOSIC 

^BSLYP 

1—1 

N2 

15,56 

10,43 

18,21 

11,8 

15,65 

O 2 

12,07 

6,92 

16,06 

7,01 

12,74 

H 2 S 

10,46 

6,36 

11,69 

7,13 

11.01 

CO 

14,10 

9,12 

15,75 

10.40 

14.10 

C 02 

13,77 

9.32 

16,26 

10.32 

14.04 

H 20 

12,62 

7,33 

15,42 

8.51 

13.11 

HCN 

13,60 

9.18 

15,41 

10.07 

14.11 

LiCl 

10,01 

6.00 

11,52 

6.76 

10.33 

LiH 

7,9 

4.39 

9.16 

5.28 

8.19 

MAE 


4.6 

2.2 

3.6 

0.4 
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TABLE II: Test set 2: Experimental values for the experimental ionization potential (IP) and the 
respective calculated negative of the HOMO energy for different methods (all values in [eV]). MAE 
is the mean average error. 


molecule 

exp. IP 

^LDA 

^FOSIC 

^ADSIC 

^BSLYP 

j-pLDA 

^SCF 

ethan 

11,52 

8.05 

14,74 

- 

9.31 

11.74 

ethen 

10,51 

6.97 

12,63 

12.1 

7.56 

11.02 

ethen — ffourene 

10,12 

6.37 

13.55 

- 

7.36 

10.25 

ethin 

11,40 

7.36 

12,94 

- 

8.09 

11.77 

benzen 

9,24 

6,55 

9,41 

9.96 

6.97 

9.55 

fulvene 

8,55 

5,56 

9,01 

8.91 

6.07 

8.60 

furan 

8,89 

5,85 

9,56 

9.85 

6.39 

9.19 

pyrazine 

9,4 

5,93 

10,67 

9.36 

7.00 

9.15 

pyridine 

9,6 

5,97 

10,17 

9.39 

7.07 

9.48 

pyrimidine 

9,73 

6,07 

10,79 

9.52 

7.11 

9.39 

thiophene 

8,87 

6,06 

9,68 

9.75 

6.56 

9.19 

dimethylsilaethene 

8,3 

5,06 

9,64 

8.6 

5.62 

8.19 

TCNQ 

7,29 

7,24 

11,80 

- 

7.56 

9.32 

naphthalene 

8,14 

5.69 

8,73 

- 

6.03 

8.16 

anthracene 

7,44 

5,19 

8,03 

- 

5.44 

7.33 

chrysene 

7,60 

5.48 

8.27 

- 

5.73 

7.45 

picene 

7,51 

5.47 

8.15* 

- 

5.71 

7.28 

MAE 


2.9 

1,4 


2.3 

0.3 


a 0.1 eV lower IP than chrysene. 

We also included some values from the literature obtained by the average density SIC ap¬ 
proach [21]. The values agree very well with our results. The average density SIC results are 
based on the BeckeSS GCA functional m which somewhat limits the direct comparison to 
our PW92 FOSIC values. However it supports the correctness of the FOSIC implementation 
and makes us optimistic to see further improvement of orbital eigenvalues when applying 
FOSIC to GGA type functionals. 
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FIG. 2: Comparison of the ionization potentials obtained from orbital eigenvalues and experimental 
values. 

Further, we found that the positions of the Fermi orbitals are transferable between similar 
molecules as previously reported [T]. It was for example very successful to use the Fermi 
orbital positions obtained for benzene as starting values for the 6-rings of the poly-aromatic 
hydrocarbons naphthalene to picene. The subsequent optimization of such an initial guess 
converged almost immediately and the resulting — ch values are in very good agreement 
with the experimental data. In Fig. [^we show such a case. For the case of naphthalene 
and anthracene one sees that the relative positions of the Fermi orbitals with respect to 
the neighbouring carbon atoms within the aromatic rings are almost identical. The Fermi 
orbitals representing single bonds are always located in the xy-plane in the middle between 
the neighbouring C atoms. On the other hand, the Fermi orbitals representing the double 
bonds are 0.3 and 1.4 A above the molecular plane. The same argument holds for the 
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Naphthalene Anthracene 

s ^ ^ s 



FIG. 3: Optimized positions of the Fermi orbitals for different structures. The red dots represent 
the final Fermi orbital positions. 

comparison of the structurally similar thiophene and furan molecules (Fig. [^bottom). Again 
the distribution of the Fermi orbitals inside the 5-ring is nearly identical. The difference 
between furan and thiophene can be identified in the additional electrons of the sulfur atom, 
that form a perfect tetrahedron and the distance of the Fermi orbitals forming the free 
electron pairs (0.90 A sulfur vs 0.78 A oxygen). The larger distance in case of thiophene is 
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FIG. 4: Position of the orbital eigenvalues of PTCDA from LDA and FOSIC method compared 
to experimental data (green) obtained by gas phase photoelectron spectroscopy (values extracted 
from [IS] ) together with a view of an isosurface plot of the respective molecular orbitals as obtained 
by the LDA and FOSIC calculations (HOMO-2 - HOMO-4 not shown for clarity). HOMO levels 
are aligned for better comparsion. 

an indication of the larger Coulomb repulsion from the additional core electrons. The out 
of plane Fermi-Orbital positions in the systems with six-membered carbon rings have also 
been reported by Pederson and Baruah in applications to Benzene [H] 

In some cases we even obtain a more realistic level ordering from the FOSIC calculation 
than from LDA. One representative example is given in Fig. |4| Here we computed the 
electronic structure of PTCDA [52], a molecule that consists of an aromatic core and oxygen 
end groups whose structure is also shown in Fig. If one compares the calculated to 
experimental values, obtained by gas phase photoemission spectroscopy, one clearly sees that 
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LDA offers not at all an appropriate description of the measured values. Hence the orbital 
energies calculated by FOSIC are in much better agreement. The orbital analysis shows 
that within LDA the orbitals delocalized over the central aromatic ring system form the 
HOMO and HOMO-5 and they agree reasonable well with experiment. The LDA HOMO-1 
to HOMO-4 orbitals however are localized at the oxygen end groups yielding LDA energies 
that are signihcant too high due the self-interaction error and therefore they appear wrongly 
in between HOMO and HOMO-5. With FOSIC the energy of these orbitals is shifted down 
below the delocalized HOMO-5 resulting in a level ordering in very good agreement with the 
experiment. The delocalized orbitals do suffer much less from the spurious self-interaction 
because the contributions to the Coulomb-potential come from a much larger distance. In 
general one has to state that in cases where delocalized states interact substantially with 
localized ones FOSIC delivers a qualitatively correct picture while standard LDA (and also 
GCA) probably does not. 

IV. SUMMARY AND OUTLOOK 

We have presented a validation of the FOSIC approach used in DFT. We have tested 
the performance of the approach on the direct estimation of the ionization potential from 
the orbital eigenvalues. We obtain ionization potentials that are comparable to the ASCF 
results and therefore in very good agreement with experiments for all molecular systems. 
We have shown that FOSIC offers a significant improvement over standard LDA and B3LYP 
calculated ionization potentials. 

The investigations suggests that FOSIC, in its current state, is capable of handling or¬ 
ganic molecules in a straightforward way. However there is still room for improvement. A 
crucial point is the improvement of the initial guess for the FO positions. One possibility 
may be using the Edmiston and Ruedenberg criteria j50], which uses the Coulomb self¬ 
interaction instead of the quadratic spread for the construction of the maximally localized 
Wannier orbitals. This is however far more demanding from an implementations point of 
view Ell and is therefore left for further studies. The observation of transferable positions 
in structurally similar molecules appears to be very promising for application of the method 
to large molecules. The current implementation is based on LSDA, therefore an extension to 
other functionals in particular to GCA ones is highly desirable. A longer term goal might be 
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to design functionals that are explicitly corrected for self-interaction error using a systematic 
nnitarily invariant FOSIC methodology discussed here. 

Acknowledgments 

S.L. wants to thank the Friedrich-Ebert-Stiftnng for the financial support. Further we 
would like to thank the ZIH Dresden for providing computational resources. M.R.P thanks 
the Texas Advanced Computing Center (TACC) and National Energy Research Scientific 
Computing (NERSC) for computer time to test the minimization methods. 


[1] M. R. Pederson, A. Ruzsinszky, and J. P. Perdew, The Journal of Chemical Physics 140, 
121103 (2014). 

[2] M. R. Pederson, J. Chem. Phys. 142, 064112 (2015). 

[3] A. J. Cohen, P. Mori-Sanchez, and W. Yang, Chem. Rev. 112, 289 (2012). 

[4] R. M. Martin, Electronic structure - basic theory and practical methods (Cambridge University 
Press, 2004). 

[5] A. D. Becke, Journal of Chemical Physics 98, 1372 (1993). 

[6] J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118 (2003). 

[7] Y. Zhao and D. Truhlar, Theoretical Chemistry Accounts 120, 215 (2008). 

[8] M. Griining, O. Gritsenko, S. Van Gisbergen, and E. Baerends, The Journal of Physical 
Chemistry A 105, 9211 (2001). 

[9] J. Kortus, C. S. Hellberg, and M. R. Pederson, Physical review letters 86, 3400 (2001). 

[10] K. Park, M. R. Pederson, and C. S. Hellberg, Physical Review B 69, 014416 (2004). 

[11] M. R. Pederson, R. A. Heaton, and C. C. Lin, The Journal of Chemical Physics 82, 2688 
(1985). 

[12] A. Dreuw, J. L. Weisman, and M. Head-Cordon, The Journal of chemical physics 119, 2943 
(2003). 

[13] T. Baruah and M. R. Pederson, Journal of Chemical Theory and Computation 5, 834 (2009). 

[14] E. Ruiz, D. R. Salahub, and A. Vela, J. Phys. Chem. 100, 12265 (1996). 

[15] S. Kiimmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). 


15 



[16] E. Ruiz, P. Alemany, S. Alvarez, and J. Cano, Journal of the American Chemical Society 
119 , 1297 (1997). 

[17] J. P. Perdew, Phys. Rev. B 23, 5048 (1981). 

[18] A. Svane and O. Gunnarsson, EPL (Europhysics Letters) 7 , 171 (1988). 

[19] B. G. Johnson, C. A. Gonzales, P. M. Gill, and J. A. Pople, Chem. Phys. Lett. 221, 100 
(1994). 

[20] S. Kliipfel, P. Kliipfel, and H. Jonsson, J. Chem. Phys. 137 , 124102 (2012). 

[21] 1. Ciofini, H. Chermette, and C. Adamo, Chem. Phys. Lett. 380, 12 (2003). 

[22] T. Tsuneda and K. Hirao, J. Chem. Phys. 140, 18A513 (2014). 

[23] T. Koopmans, Physica 1 , 104 (1934). 

[24] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982). 

[25] R. A. Heaton, J. G. Harrison, and C. C. Lin, Physical Review B 31 , 1077 (1985). 

[26] 1. Dabo et a/., Phys. Rev. B 82, 115121 (2010). 

[27] 1 . Dabo et al.^ Physical Chemistry Chemical Physics 15 , 685 (2013). 

[28] J. Perdew and A. Zunger, Physical Review B 23 , 5048 (1981). 

[29] R. A. Heaton, J. G. Harrison, and C. C. Lin, Physical Review B 28 , 5992 (1983). 

[30] R. A. Heaton, M. R. Pederson, and C. C. Lin, The Journal of chemical physics 86, 258 (1987). 

[31] M. R. Pederson and C. C. Lin, The Journal of chemical physics 88, 1807 (1988). 

[32] W. L. Luken and D. N. Beratan, Theoretica chimica acta 61, 265 (1982). 

[33] W. L. Luken and J. C. Culberson, Theoretica chimica acta 66, 279 (1984). 

[34] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997). 

[35] M. Pederson and J. Perdew, Psi-K Scientific Highlight of the Month (2012). 

[36] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. Lett. 94 , 026405 (2005). 

[37] E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, and P. Gumbsch, Phys. Rev. Lett. 97 , 170201 
(2006). 

[38] K. Jackson and M. Pederson, Physical Review B 42 , 3276 (1990). 

[39] M. Pederson, D. Porezag, J. Kortus, and D. Patton, physica status solidi (b) 217 , 197 (2000). 

[40] D. Porezag and M. R. Pederson, Physical Review A 60, 2840 (1999). 

[41] M. R. Pederson and K. A. Jackson, Phys. Rev. B 41, 7453 (1990). 

[42] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 

[43] A. Schafer, H. Horn, and R. Ahlrichs, The Journal of Chemical Physics 97 , 2571 (1992). 


16 



[44] F. Weigend and R. Ahlrichs, Physical Chemistry Chemical Physics 7 , 3297 (2005). 

[45] F. Neese, International Journal of Quantum Chemistry 83, 104 (2001). 

[46] F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7 , 3297 (2005). 

[47] A. D. Becke, Phys. Rev. A 38, 3098 (1988). 

[48] M. R. Pederson and T. Baruah, Advances in Atomic, Molecular and Optical Physics 64 , 207 
(2015). 

[49] N. Dori et a/., Phys. Rev. B 73 , 195208 (2006). 

[50] C. Edmiston and K. Ruedenberg, Rev. Mod. Phys. 35, 457 (1963). 

[51] J. E. Subotnik, Y. Shao, W. Liang, and M. Head-Cordon, J. Chem. Phys. 121, 9220 (2004). 

[52] 3,4,9,10-Perylenetetracarboxylic dianhydride 


17 



